home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
c
/
cuj9205.zip
/
1005032A
< prev
next >
Wrap
Text File
|
1992-06-02
|
2KB
|
52 lines
#include <math.h>
/* Adaptive quadrature routine.
Called with:
lft: left endpoint
rt: right endpoint
f: pointer to function to integrate
err: pointer to allowable error
Returns
return value evaluated integral
err pointer to upper limit of actual error
*/
double adaptive (double lft, double rt, double (*f)(double x), double *err)
{
double h; /* Width of this interval. */
double p1, p2, p3, p4, p5; /* Value of function at the 5 points. */
double error; /* Actual error on this interval. */
double err1, err2; /* Errors of the sub-intervals. */
double int3, /* Value of the integral with Simpson */
int5; /* 3- and 5-point approximations. */
h = rt - lft; /* Get the width of the interval. */
/* Get the 2 approximations to the integral. */
p1 = (*f)(lft);
p2 = (*f)(lft + h/4.0);
p3 = (*f)(lft + h/2.0);
p4 = (*f)(rt - h/4.0);
p5 = (*f)(rt);
int3 = h * (p1 + 4.0 * p3 + p5) / 6.0;
int5 = h * (p1 + 4.0 * p2 + 2.0 * p3 + 4.0 * p4 + p5) / 12.0;
error = fabs (int3 - int5) / 15.0; /* Find the actual error. */
if (error < *err) /* If within tolerance, */
{
*err = error; /* copy the actual error to caller, */
return (int5); /* return the more accurate approx. */
}
else /* Otherwise, */
{
err1 = err2 = *err / 2.0; /* cut the allowable error in half, */
int5 = adaptive (lft, lft + h/2.0, f, &err1); /* integrate the two */
int5 += adaptive (lft + h/2.0, rt, f, &err2); /* half-regions, */
*err = err1 + err2; /* sum the real errors, */
return (int5); /* and return the integral. */
}
}